home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / log1p.c < prev    next >
C/C++ Source or Header  |  1988-07-11  |  5KB  |  163 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms are permitted
  6.  * provided that this notice is preserved and that due credit is given
  7.  * to the University of California at Berkeley. The name of the University
  8.  * may not be used to endorse or promote products derived from this
  9.  * software without specific prior written permission. This software
  10.  * is provided ``as is'' without express or implied warranty.
  11.  *
  12.  * All recipients should regard themselves as participants in an ongoing
  13.  * research project and hence should feel obligated to report their
  14.  * experiences (good or bad) with these elementary function codes, using
  15.  * the sendbug(8) program, to the authors.
  16.  */
  17.  
  18. #ifndef lint
  19. static char sccsid[] = "@(#)log1p.c    5.2 (Berkeley) 4/29/88";
  20. #endif /* not lint */
  21.  
  22. /* LOG1P(x) 
  23.  * RETURN THE LOGARITHM OF 1+x
  24.  * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS)
  25.  * CODED IN C BY K.C. NG, 1/19/85; 
  26.  * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85.
  27.  * 
  28.  * Required system supported functions:
  29.  *    scalb(x,n) 
  30.  *    copysign(x,y)
  31.  *    logb(x)    
  32.  *    finite(x)
  33.  *
  34.  * Required kernel function:
  35.  *    log__L(z)
  36.  *
  37.  * Method :
  38.  *    1. Argument Reduction: find k and f such that 
  39.  *            1+x  = 2^k * (1+f), 
  40.  *       where  sqrt(2)/2 < 1+f < sqrt(2) .
  41.  *
  42.  *    2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
  43.  *         = 2s + 2/3 s**3 + 2/5 s**5 + .....,
  44.  *       log(1+f) is computed by
  45.  *
  46.  *                 log(1+f) = 2s + s*log__L(s*s)
  47.  *       where
  48.  *        log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...)))
  49.  *
  50.  *       See log__L() for the values of the coefficients.
  51.  *
  52.  *    3. Finally,  log(1+x) = k*ln2 + log(1+f).  
  53.  *
  54.  *    Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers
  55.  *           n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last 
  56.  *           20 bits (for VAX D format), or the last 21 bits ( for IEEE 
  57.  *           double) is 0. This ensures n*ln2hi is exactly representable.
  58.  *        2. In step 1, f may not be representable. A correction term c
  59.  *            for f is computed. It follows that the correction term for
  60.  *           f - t (the leading term of log(1+f) in step 2) is c-c*x. We
  61.  *           add this correction term to n*ln2lo to attenuate the error.
  62.  *
  63.  *
  64.  * Special cases:
  65.  *    log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal;
  66.  *    log1p(INF) is +INF; log1p(-1) is -INF with signal;
  67.  *    only log1p(0)=0 is exact for finite argument.
  68.  *
  69.  * Accuracy:
  70.  *    log1p(x) returns the exact log(1+x) nearly rounded. In a test run 
  71.  *    with 1,536,000 random arguments on a VAX, the maximum observed
  72.  *    error was .846 ulps (units in the last place).
  73.  *
  74.  * Constants:
  75.  * The hexadecimal values are the intended ones for the following constants.
  76.  * The decimal values may be used, provided that the compiler will convert
  77.  * from decimal to binary accurately enough to produce the hexadecimal values
  78.  * shown.
  79.  */
  80.  
  81. #if defined(vax)||defined(tahoe)    /* VAX D format */
  82. #include <errno.h>
  83. #ifdef vax
  84. #define _0x(A,B)    0x/**/A/**/B
  85. #else    /* vax */
  86. #define _0x(A,B)    0x/**/B/**/A
  87. #endif    /* vax */
  88. /* static double */
  89. /* ln2hi  =  6.9314718055829871446E-1    , Hex  2^  0   *  .B17217F7D00000 */
  90. /* ln2lo  =  1.6465949582897081279E-12   , Hex  2^-39   *  .E7BCD5E4F1D9CC */
  91. /* sqrt2  =  1.4142135623730950622E0     ; Hex  2^  1   *  .B504F333F9DE65 */
  92. static long     ln2hix[] = { _0x(7217,4031), _0x(0000,f7d0)};
  93. static long     ln2lox[] = { _0x(bcd5,2ce7), _0x(d9cc,e4f1)};
  94. static long     sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)};
  95. #define    ln2hi    (*(double*)ln2hix)
  96. #define    ln2lo    (*(double*)ln2lox)
  97. #define    sqrt2    (*(double*)sqrt2x)
  98. #else    /* defined(vax)||defined(tahoe)    */
  99. static double
  100. ln2hi  =  6.9314718036912381649E-1    , /*Hex  2^ -1   *  1.62E42FEE00000 */
  101. ln2lo  =  1.9082149292705877000E-10   , /*Hex  2^-33   *  1.A39EF35793C76 */
  102. sqrt2  =  1.4142135623730951455E0     ; /*Hex  2^  0   *  1.6A09E667F3BCD */
  103. #endif    /* defined(vax)||defined(tahoe)    */
  104.  
  105. double log1p(x)
  106. double x;
  107. {
  108.     static double zero=0.0, negone= -1.0, one=1.0, 
  109.               half=1.0/2.0, small=1.0E-20;   /* 1+small == 1 */
  110.     double logb(),copysign(),scalb(),log__L(),z,s,t,c;
  111.     int k,finite();
  112.  
  113. #if !defined(vax)&&!defined(tahoe)
  114.     if(x!=x) return(x);    /* x is NaN */
  115. #endif    /* !defined(vax)&&!defined(tahoe) */
  116.  
  117.     if(finite(x)) {
  118.        if( x > negone ) {
  119.  
  120.        /* argument reduction */
  121.           if(copysign(x,one)<small) return(x);
  122.           k=logb(one+x); z=scalb(x,-k); t=scalb(one,-k);
  123.           if(z+t >= sqrt2 ) 
  124.           { k += 1 ; z *= half; t *= half; }
  125.           t += negone; x = z + t;
  126.           c = (t-x)+z ;        /* correction term for x */
  127.  
  128.         /* compute log(1+x)  */
  129.               s = x/(2+x); t = x*x*half;
  130.           c += (k*ln2lo-c*x);
  131.           z = c+s*(t+log__L(s*s));
  132.           x += (z - t) ;
  133.  
  134.           return(k*ln2hi+x);
  135.        }
  136.     /* end of if (x > negone) */
  137.  
  138.         else {
  139. #if defined(vax)||defined(tahoe)
  140.         extern double infnan();
  141.         if ( x == negone )
  142.             return (infnan(-ERANGE));    /* -INF */
  143.         else
  144.             return (infnan(EDOM));    /* NaN */
  145. #else    /* defined(vax)||defined(tahoe) */
  146.         /* x = -1, return -INF with signal */
  147.         if ( x == negone ) return( negone/zero );
  148.  
  149.         /* negative argument for log, return NaN with signal */
  150.             else return ( zero / zero );
  151. #endif    /* defined(vax)||defined(tahoe) */
  152.         }
  153.     }
  154.     /* end of if (finite(x)) */
  155.  
  156.     /* log(-INF) is NaN */
  157.     else if(x<0) 
  158.          return(zero/zero);
  159.  
  160.     /* log(+INF) is INF */
  161.     else return(x);      
  162. }
  163.